Recent advances in computer industry made possible a wide variety of biological molecular modeling applications with explicit treatment of atomic sites on the molecules of interest. Computing pKa's of proteins is one of the important tasks crucial to allow construction of the protonation states desired by varying the protein and solvent environment and thus to facilitate the contemporary drug design. This problem has been approached by many researchers via modeling proteins (or part of them) at the explicit all-atom level and representing the solvent with the dielectric continuum. However there is still a significant discrepancy between the theoretical and experimental results, as well as a controversy in the issue of selecting the most adequate computational method for this task. Two major problems are (i) the problem of computing the mani-body energy contribution and (ii) choosing the correct geometries for the solvated molecules. Commonly assigned to a protein "effective dielectric constant" does not allow results much better than those given by the trivial null-model. On the other hand, using crystallographically obtained molecular geometries in solution for both protonated and un-protonated states is one of the most important sources on the discrepancy between the experimental and theoretical results as well. We will conduct geometry optimizations of model systems and real proteins in a dielectric continuum solvent (water). Geometries of both protonated and un-protonated states will be optimized. We plan to abandon the "effective dielectric constant" of the proteins and to employ the induced point dipole polarization interactions to represent the many-body energy. Energetic results of the geometry optimizations will be used to compute pKa's. We hope that this will let us to obtain pKa's of proteins in better agreement with the experimental data and will shade light on mechanisms on the microscopic processes in biochemical systems.